Sistemas y Señales Biomédicos

Ingeniería Biomédica

Ph.D. Pablo Eduardo Caicedo Rodríguez

2025-05-05

Sistemas y Señales Biomedicos - SYSB

Digital Filters

Why the Z-Transform?

  • The Fourier Transform assumes signals are stable and well-behaved
  • But some biosignals or systems may not be absolutely summable
  • The Z-Transform generalizes the Fourier Transform
  • Useful for analyzing discrete-time systems, especially when stability and causality matter

Definition

Let \(x[n]\) be a discrete-time signal.

The Z-Transform is defined as:

\[X(z) = \sum_{n=-\infty}^{\infty} x[n] z^{-n}\]

Where: - \(z \in \mathbb{C}\) is a complex variable - \(z = re^{j\omega}\)

Region of Convergence (ROC)

  • The Z-Transform converges only for certain values of \(z\)
  • The set of \(z\) for which the series converges is the ROC
  • ROC is critical for system stability and causality

Causal Signals
ROC is outside outermost pole

Anti-Causal Signals
ROC is inside innermost pole

Z-Plane Representation

  • Poles: values of \(z\) where \(X(z) \to \infty\)
  • Zeros: values where \(X(z) = 0\)
  • Visualization of poles and zeros helps in understanding system behavior

\[H(z) = 1.00 \cdot \frac{(z - 0.50)}{(z - 0.90)}\]

(-1.5, 1.5)
(-1.5, 1.5)

Relationship with Fourier Transform

If the ROC includes the unit circle, \(|z| = 1\), then:

\[X(e^{j\omega}) = \sum_{n=-\infty}^{\infty} x[n] e^{-j\omega n}\]

So the Fourier Transform is a special case of the Z-Transform.

Example Signal

Let:

\[x[n] = a^n u[n]\]

Where: - \(a \in \mathbb{R}\) - \(u[n]\) is the unit step function (0 for \(n<0\), 1 for \(n\geq 0\))

Step 1: Z-Transform

Apply the definition:

\[X(z) = \sum_{n=0}^{\infty} a^n z^{-n} = \sum_{n=0}^{\infty} (az^{-1})^n\]

This is a geometric series:

\[X(z) = \frac{1}{1 - az^{-1}} = \frac{z}{z - a}\]

What is a Geometric Series?

A geometric series is a sum of terms where each term is multiplied by the same constant:

\[S = \sum_{n=0}^{\infty} r^n = 1 + r + r^2 + r^3 + \cdots\]

The value of \(r\) determines whether this sum converges (has a finite limit) or diverges.

Goal

Understand why this series:

\[\sum_{n=0}^{\infty} r^n\]

converges if and only if:

\[\boxed{|r| < 1}\]

Partial Sums

Let’s consider the sum up to the \(N\)-th term:

\[S_N = \sum_{n=0}^{N} r^n = 1 + r + r^2 + \cdots + r^N\]

This has a closed-form expression:

\[S_N = \frac{1 - r^{N+1}}{1 - r} \quad \text{for } r \neq 1\]

Taking the Limit

To find the sum of the infinite series, take the limit as \(N \to \infty\):

\[S = \lim_{N \to \infty} S_N = \lim_{N \to \infty} \frac{1 - r^{N+1}}{1 - r}\]

Case 1: \(|r| < 1\)

If \(|r| < 1\), then:

\[r^{N+1} \to 0 \quad \text{as } N \to \infty\]

So the sum becomes:

\[S = \frac{1}{1 - r}\]

✅ The geometric series converges.

Case 2: \(|r| \geq 1\)

  • If \(r = 1\), then \(S_N = N + 1 \to \infty\)
  • If \(r = -1\), the sum oscillates: \(1 - 1 + 1 - 1 + \cdots\)
  • If \(|r| > 1\), then \(r^{N+1} \to \infty\)

❌ In all cases: the series diverges

Intuition

  • When \(|r| < 1\), each term \(r^n\) gets smaller and smaller
  • Their total sum settles to a finite number
  • When \(|r| \geq 1\), the terms don’t vanish — the sum keeps growing or oscillating

Example: \(r = 0.5\)

\[S = 1 + 0.5 + 0.25 + 0.125 + \cdots = \frac{1}{1 - 0.5} = 2\]

Every term adds less. The sum “flattens out.”

Why This Matters

The Z-Transform often gives us geometric series like:

\[\sum_{n=0}^{\infty} (az^{-1})^n\]

This converges only if:

\[|az^{-1}| < 1 \Rightarrow |z| > |a|\]

So, understanding convergence of geometric series = understanding ROC in Z-transforms.

Summary

Condition Behavior Result
\(|r| < 1\) Terms shrink Series converges
\(|r| \geq 1\) Terms grow or oscillate Diverges

\[\sum_{n=0}^{\infty} r^n = \frac{1}{1 - r} \quad \text{if } |r| < 1\]

Step 2: Region of Convergence (ROC)

For convergence of the geometric series:

\[|az^{-1}| < 1 \Rightarrow |z| > |a|\]

Therefore, the ROC is:

\[\boxed{|z| > |a|}\]

Interpretation

  • Causal signal (defined for \(n \geq 0\))
  • ROC is outside the outermost pole
  • Stable system only if ROC includes the unit circle \(|z| = 1\)

Example: \(a = 0.5\)

\[x[n] = (0.5)^n u[n]\]

Z-Transform:

\[X(z) = \sum_{n=0}^{\infty} (0.5)^n z^{-n} = \frac{z}{z - 0.5}\]

Region of Convergence:

\[\boxed{|z| > 0.5}\]

Summary

Signal Z-Transform ROC
\(x[n] = a^n u[n]\) \(\frac{z}{z - a}\) \(|z| > |a|\)
Example: \(a = 0.5\) \(\frac{z}{z - 0.5}\) \(|z| > 0.5\)

Properties of the Z-Transform

  • Linearity: \(a x[n] + b y[n] \to aX(z) + bY(z)\)
  • Time shifting: \(x[n - k] \to z^{-k} X(z)\)
  • Scaling in the z-domain: \(a^n x[n] \to X(z/a)\)
  • Convolution: \(x[n] * h[n] \to X(z)H(z)\)

Example

Let \(x[n] = a^n u[n]\), where \(|a| < 1\)

\[X(z) = \sum_{n=0}^{\infty} a^n z^{-n} = \frac{1}{1 - az^{-1}}, \quad \text{ROC: } |z| > |a|\]

Difference Equations in DSP

A difference equation relates input and output values at different time steps.

\[y[n] - a_1 y[n-1] - a_2 y[n-2] = b_0 x[n] + b_1 x[n-1]\]

Common in: - Digital filters (FIR, IIR) - Signal models in ECG, EEG analysis - Implementation in real-time biosignal systems

Z-Transform of Time-Shifted Terms

The Z-Transform turns time shifts into powers of \(z^{-1}\):

Time Domain Z-Domain
\(x[n]\) \(X(z)\)
\(x[n-k]\) \(z^{-k} X(z)\)
\(y[n-k]\) \(z^{-k} Y(z)\)

Step 1: Apply Z-Transform

Given:

\[y[n] - a_1 y[n-1] - a_2 y[n-2] = b_0 x[n] + b_1 x[n-1]\]

Apply \(\mathcal{Z} \{ \cdot \}\):

\[Y(z) - a_1 z^{-1} Y(z) - a_2 z^{-2} Y(z) = b_0 X(z) + b_1 z^{-1} X(z)\]

Step 2: Factor and Solve for \(H(z)\)

Group:

\[Y(z)(1 - a_1 z^{-1} - a_2 z^{-2}) = X(z)(b_0 + b_1 z^{-1})\]

Divide both sides:

\[H(z) = \frac{Y(z)}{X(z)} = \frac{b_0 + b_1 z^{-1}}{1 - a_1 z^{-1} - a_2 z^{-2}}\]

Example

Given:

\[y[n] - 0.9 y[n-1] = x[n] - 0.5 x[n-1]\]

Z-Transform:

\[Y(z)(1 - 0.9 z^{-1}) = X(z)(1 - 0.5 z^{-1})\]

Transfer Function:

\[H(z) = \frac{1 - 0.5 z^{-1}}{1 - 0.9 z^{-1}}\]

Poles and Zeros

Let’s analyze \(H(z)\):

  • Zeros: Roots of the numerator \(\Rightarrow z = 0.5\)
  • Poles: Roots of the denominator \(\Rightarrow z = 0.9\)

Pole-Zero Plot
Visualizes system behavior
Check for: - Stability (poles inside unit circle) - Frequency shaping

Practice

Convert this equation:

\[y[n] = 0.6 y[n-1] + x[n] + x[n-1]\]

Find: - \(H(z)\) - Poles and zeros - Plot them in the Z-plane

Application in Biosignal Processing

  • Analysis of digital filters for ECG, EEG, etc.
  • Design of stable and causal filtering systems
  • Useful in difference equation modeling of biosignals

Summary

  • Z-Transform is a powerful tool for analyzing discrete systems
  • Provides insight into stability, causality, and system behavior
  • A generalization of the Fourier Transform
  • Crucial in digital signal processing of biosignals
  • Z-Transform converts difference equations into algebraic expressions
  • Transfer function \(H(z)\) tells us how the system responds to inputs
  • Key for digital filter design in biosignal processing

Next Steps

  • Practice Z-Transform computations
  • Pole-zero plotting exercises
  • Application to real biosignal filtering problems

Filter Design

Notch Filter

A common issue in biosignal processing is removing power‐line interference (50/60 Hz) from, for example, EEG or ECG signals. A simple digital filter to eliminate 60 Hz interference (assuming a sampling frequency \(f_s = 5000\) Hz) is to place complex‐conjugate zeros at

\[ z = e^{\pm j 2\pi\frac{60}{5000}}. \]

The resulting transfer function can be written as

\[ H(z) = 1 \;-\; 2\cos\!\Bigl(2\pi\frac{60}{5000}\Bigr)\,z^{-1} \;+\; z^{-2}. \]

This \(H(z)\) has zeros at \(e^{\pm j2\pi(60/5000)}\) that precisely cancel the 60 Hz component, thereby implementing a notch filter. Moreover, it is a second‐order FIR filter with symmetric coefficients, which grants it linear phase (important to avoid waveform distortion).

Filter Design

def design_filter(zeros=None, poles=None, gain=1.0):
    """
    Diseña un filtro digital a partir de ceros y/o polos y una ganancia.

    Parámetros:
    - zeros: lista de ceros (raíces del numerador), o None para no incluir
    - poles: lista de polos (raíces del denominador), o None para no incluir
    - gain: ganancia escalar del filtro

    Devuelve:
    - b: coeficientes del numerador
    - a: coeficientes del denominador
    """
    # Si no se pasan ceros, asumimos un FIR trivial (b = [gain])
    if zeros:
        b = gain * np.poly(zeros)
    else:
        b = np.array([gain], dtype=float)

    # Si no se pasan polos, asumimos sistema FIR (a = [1])
    if poles:
        a = np.poly(poles)
    else:
        a = np.array([1.0], dtype=float)

    return b, a

Filter Design

def gaussian(x, mu, sigma, A):
    """
    Genera una función gaussiana.

    Parámetros:
    - x: array de tiempos
    - mu: posición central de la gaussiana
    - sigma: desviación estándar
    - A: amplitud
    """
    return A * np.exp(-((x - mu) ** 2) / (2 * sigma**2))

Filter Design

def simulate_ecg(duration=10.0, fs=500, heart_rate=60):
    """
    Simula un ECG sintético basado en la superposición de ondas gaussianas.

    Parámetros:
    - duration: duración de la señal en segundos
    - fs: frecuencia de muestreo en Hz
    - heart_rate: frecuencia cardiaca en latidos por minuto

    Devuelve:
    - t: vector de tiempos
    - ecg: señal simulada de ECG en mV
    """
    dt = 1 / fs
    t = np.arange(0, duration, dt)
    rr = 60 / heart_rate  # intervalo RR en segundos

    # Inicializar señal
    ecg = np.zeros_like(t)

    # Parámetros de las ondas (posiciones relativas en segundos)
    # P wave
    p_amp, p_dur, p_delay = 0.25, 0.09, 0.16
    # Q wave
    q_amp, q_dur, q_delay = -0.05, 0.066, 0.166
    # R wave
    r_amp, r_dur, r_delay = 1.6, 0.1, 0.166
    # S wave
    s_amp, s_dur, s_delay = -0.25, 0.066, 0.19
    # T wave
    t_amp, t_dur, t_delay = 0.35, 0.142, 0.36

    # Generar cada latido
    for beat_start in np.arange(0, duration, rr):
        mask = (t >= beat_start) & (t < beat_start + rr)
        tb = t[mask] - beat_start
        ecg[mask] += gaussian(tb, p_delay, p_dur / 2, p_amp)
        ecg[mask] += gaussian(tb, q_delay, q_dur / 2, q_amp)
        ecg[mask] += gaussian(tb, r_delay, r_dur / 2, r_amp)
        ecg[mask] += gaussian(tb, s_delay, s_dur / 2, s_amp)
        ecg[mask] += gaussian(tb, t_delay, t_dur / 2, t_amp)

    return t, ecg+0.1*np.cos(2*np.pi*60*t)

Filter Design

# Parámetros de simulación
DURATION = 10    # segundos
FS = 500         # Hz
HR = 70          # latidos por minuto

# Generar señal
t, ecg_signal = simulate_ecg(duration=DURATION, fs=FS, heart_rate=HR)

# Graficar resultado
plt.figure(figsize=(12, 4))
plt.plot(t, ecg_signal, linewidth=1)
plt.title(f'Señal de ECG sintética ({HR} bpm)')
plt.xlabel('Tiempo (s)')
plt.ylabel('Amplitud (mV)')
plt.grid(True)
plt.tight_layout()
plt.show()

Filter Design

n = len(ecg_signal)
yf = np.fft.fft(ecg_signal)
xf = np.fft.fftfreq(n, 1/fs)[:n//2]
plt.plot(xf, 2.0/n * np.abs(yf[0:n//2]))
plt.grid()
plt.show()

Filter Design

fc = 60
b,a =design_filter(zeros=[np.exp(1j*2*np.pi*fc/fs), np.exp(-1j*2*np.pi*fc/fs)])
print(a)
print(b)

w, h = sig.freqz(a, b, worN=8000)
plt.plot(w/np.pi*fs/2, 20*np.log10(abs(h)))
plt.grid()
plt.xlabel('Frequency (Hz)')

Filter Design

FIR filters

FIR filters (Finite Impulse Response) are widely used in biomedical processing because they can be designed to have linear phase response, avoiding phase distortion in the filtered signal (which is useful for preserving the morphology of ECG waves, for example).

Filter Design

Filter Design Process

  1. Defining the desired ideal frequency response \(H_d(e^{j\omega})\).

  2. Obtaining the ideal impulse response \(h_d[n]\) as the inverse Fourier transform of \(H_d\).

  3. Truncating \(h_d[n]\) (which is usually infinite or very long) with a window function \(w[n]\) to obtain a realizable FIR filter

    \[ h[n] = h_d[n]\,w[n]. \]

Filter Design

You obtain \(h_d[n]\) as the inverse discrete–time Fourier transform of your ideal frequency response \(H_d(e^{j\omega})\). For a low-pass filter with cutoff \(\omega_c\),

\[ H_d(e^{j\omega}) = \begin{cases} 1, & |\omega|\le\omega_c,\\ 0, & \text{otherwise}. \end{cases} \]

By definition of the inverse DTFT,

\[ h_d[n] \;=\; \frac{1}{2\pi}\int_{-\pi}^{\pi} H_d(e^{j\omega})\,e^{j\omega n}\,d\omega \;=\;\frac{1}{2\pi}\int_{-\omega_c}^{\omega_c} e^{j\omega n}\,d\omega. \]

Filter Design

Carry out the integral in two cases:

  1. For \(n\neq 0\):

    \[ h_d[n] =\frac{1}{2\pi}\,\frac{e^{j\omega n}}{j\,n}\Biggr|_{-\omega_c}^{\omega_c} =\frac{1}{2\pi}\,\frac{e^{j\omega_c n}-e^{-j\omega_c n}}{j\,n} =\frac{\sin(\omega_c\,n)}{\pi\,n}. \]

  2. For \(n=0\):

    \[ h_d[0] =\frac{1}{2\pi}\int_{-\omega_c}^{\omega_c} 1\,d\omega =\frac{2\,\omega_c}{2\pi} =\frac{\omega_c}{\pi}. \]

Putting both together,

\[ h_d[n] =\begin{cases} \dfrac{\sin(\omega_c\,n)}{\pi\,n}, & n\neq 0,\\[1em] \dfrac{\omega_c}{\pi}, & n=0. \end{cases} \]

Relating to sampling frequency

If your cutoff is specified in Hz, \(f_c\), and sampling rate is \(f_s\), then

\[ \omega_c = 2\pi\,\frac{f_c}{f_s}, \]

so you can write

\[ h_d[n] =\frac{\sin\!\bigl(2\pi \frac{f_c}{f_s}\,n\bigr)}{\pi\,n} =\;2\;\frac{f_c}{f_s}\;\frac{\sin\!\bigl(2\pi \frac{f_c}{f_s}\,n\bigr)}{2\pi \frac{f_c}{f_s}\,n} =2\frac{f_c}{f_s}\,\mathrm{sinc}\!\Bigl(2\frac{f_c}{f_s}\,n\Bigr), \]

where we define the normalized sinc as \(\mathrm{sinc}(x)=\frac{\sin(\pi x)}{\pi x}\).

Key takeaway

  • \(h_d[n]\) is exactly the inverse‐DTFT of the ideal (“brick‐wall”) frequency specification.
  • It produces a sinc-shaped impulse response of infinite length.
  • Truncation (with a window) makes it realizable as a finite-length FIR.

Filter Design

The typical characteristics of classic windows are:

  • Rectangular: narrowest main lobe (width ≈ 4π/M rad) but highest sidelobes (first sidelobe ≈ –13 dB, stop-band attenuation ≈ 21 dB). It gives the steepest transition for a given filter order, at the expense of poorer stop-band rejection.

  • Bartlett (triangular): somewhat wider main lobe (≈ 8π/M), sidelobes ≈ –25 dB.

  • Hann: main lobe ≈ 8π/M, sidelobes ≈ –31 dB (better rejection than rectangular but smoother transitions).

  • Hamming: main lobe ≈ 8π/M, sidelobes ≈ –41 dB (minimum stop-band attenuation ≈ 53 dB). Very popular for its good compromise between transition width and stop-band attenuation.

  • Blackman: wider main lobe (≈ 12π/M) but very low sidelobes (≈ –57 dB, attenuation ≈ 74 dB).

  • Kaiser: allows selection of a parameter β to control sidelobe attenuation, offering flexibility. Approximately, to achieve A dB of attenuation,

    \[ \beta \approx 0.1102\,(A - 8.7)\quad(\text{for }A>50), \]

    and the normalized transition width Δω relates to the order M and β by

    \[ M \approx \frac{A - 8}{2.285\,\Delta\omega} \]

These formulas stem from Kaiser’s approximations and help in sizing the filter.

Windows Forms

(-100.0, 5.0)

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal.windows import (
    boxcar,       # Rectangular
    bartlett,    # Triangular
    hann,        # Hann
    hamming,     # Hamming
    blackman,    # Blackman
    kaiser       # Kaiser
)

# Parameters
M = 64               # window length
beta = 8.6           # Kaiser parameter for moderate sidelobe attenuation
nfft = 512           # FFT length for frequency response
fs = 1.0             # normalized sampling rate

# Generate windows
windows = {
    'Rectangular': boxcar(M),
    'Bartlett':    bartlett(M),
    'Hann':        hann(M),
    'Hamming':     hamming(M),
    'Blackman':    blackman(M),
    f'Kaiser (β={beta})': kaiser(M, beta)
}

# Time-domain plot
plt.figure(figsize=(8, 4))
for name, w in windows.items():
    plt.plot(np.arange(M), w, label=name)
plt.title('Window Functions — Time Domain')
plt.xlabel('Sample index n')
plt.ylabel('Amplitude')
plt.legend(loc='upper right')
plt.grid(True)
plt.tight_layout()

# Frequency-domain plot
plt.figure(figsize=(8, 4))
for name, w in windows.items():
    # Compute normalized frequency response
    W = np.fft.fft(w, nfft)
    W_mag = 20 * np.log10(np.abs(W) / np.max(np.abs(W)))
    freqs = np.fft.fftfreq(nfft, d=1/fs)
    # Only plot 0 ≤ f ≤ 0.5 (normalized Nyquist)
    idx = np.logical_and(freqs >= 0, freqs <= 0.5)
    plt.plot(freqs[idx], W_mag[idx], label=name)
plt.title('Window Functions — Frequency Response')
plt.xlabel('Normalized Frequency (cycles/sample)')
plt.ylabel('Magnitude (dB)')
plt.ylim(-100, 5)
plt.legend(loc='upper right')
plt.grid(True)
plt.tight_layout()

plt.show()

Filter Design

N = 101
fc = 30

# 2. Compute the ideal impulse response h_d[n] of a low-pass filter
n = np.arange(N)
M = (N - 1) / 2
# Using the normalized sinc: sinc(x) = sin(pi*x)/(pi*x)
h_d = (2 * fc / FS) * np.sinc(2 * fc * (n - M) / FS)

# 3. Choose a window w[n] (here: Hamming) and truncate
w = np.hamming(N)
h = h_d * w

# 4. Normalize to ensure unity gain at DC
h /= np.sum(h)

# 6. Filter it (only allowed library call)
y = sig.lfilter(h, 1.0, ecg_signal)

# 7. Plot input vs. output
plt.figure(figsize=(8, 4))
plt.plot(t, ecg_signal, label='Original signal')
plt.plot(t, y, label='Filtered signal', linewidth=2)
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.title('Low-pass FIR (window method) – Cutoff 100 Hz')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()